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We present the^r^f numerical-relativity simulation of a compact-object binary whose gravitational waveform 
is long enough to cover the entire frequency band of advanced gravitational-wave detectors, such as LIGO, 

Virgo and KAGRA, for mass ratio 7 and total mass as low as 45.5M0. We find that effective-one-body models, 
either uncalibrated or calibrated against substantially shorter numerical-relativity waveforms at smaller mass 
ratios, reproduce our new waveform remarkably well, with a negligible loss in detection rate due to model¬ 
ing error. In contrast, post-Newtonian inspiral waveforms and existing calibrated phenomenological inspiral- 
merger-ringdown waveforms display greater disagreement with our new simulation. The disagreement varies 
substantially depending on the specific post-Newtonian approximant used. 

PACS numbers: 04.25.D-, 04.25.dg, 04.30.-w, 04.30.Db 


Introduction. The upgraded ground-based interferometric 
gravitational-wave (GW) detectors LIGO [1,2] and Virgo [3] 
will begin scientific observations in mid 2015, and are ex¬ 
pected to reach design sensitivity by 2019 [4]. Furthermore, 
a new Japanese detector, KAGRA [5], is under construction. 
Direct detection of GWs by the end of this decade is therefore 
very likely. The most promising GW sources are compact- 
object binaries, wherein each partner is either a stellar-mass 
black hole (BH) or a neutron star (NS) [6]. The detection 
of GWs from compact-object binaries, as well as the deter¬ 
mination of source properties from detected GW signals, re¬ 
lies on the accurate knowledge of the expected gravitational 
waveforms via matched-filtering [7] and Markov chain Monte 
Carlo techniques. 

The need for accurate waveforms has motivated intense 
research. Early waveform models, based on the post- 
Newtonian (PN) formalism [8], were limited to the early 
inspiral. The effective-one-body (EOB) formalism [9, 10] 
extended waveform models to the late inspiral, merger and 
ringdown. Since 2005 research has greatly benefited from 
numerical-relativity (NR) simulations [11-13]*. Current 
inspiral-merger-ringdown (IMR) waveform models [17-20] 
combine information from analytical-relativity (AR) calcula¬ 
tions (best suited for the inspiral, when comparable-mass bi¬ 
naries have characteristic velocities smaller than the speed of 
light) and direct NR simulations (the best means to explore 
the late inspiral and the merger). 


However, there is a gap between the portion of the binary 
evolution that is described by analytical methods, and the por¬ 
tion that is accessible by NR. For example, waveforms com¬ 
puted at the currently-known PN order become unreliable pos¬ 
sibly hundreds of orbits before merger for unequal-mass bi¬ 
naries [21, 22], and even earlier when one of the objects is 
spinning [23] ^ yet NR simulations have been able to cover 
only tens of orbits [24-26] until now. This gap has emerged 
as perhaps the most important source of uncertainty in present 
IMR waveform models. It is possible to construct IMR mod¬ 
els by extending analytical waveforms across the gap [17-20], 
in some cases obtaining IMR models that are faithful to longer 
numerical waveforms when extrapolated beyond their limited 
range of calibration [27]. However, so far these procedures 
have been tested using NR simulations with only 30 orbits, 
too few to close the gap. 

The time duration T of an inspiral waveform starting at 
initial GW frequency /jni scales as T I' where 

V = mimi/M^ is the symmetric mass ratio of the binary with 
component masses m\ 2 and total mass M — ni\ A-m 2 . There¬ 
fore, reducing /ini by a factor of 2 increases T sevenfold. 
Halving the symmetric mass ratio V (e.g., from mi jm 2 = 2 
to mi/m 2 =7) doubles T. Increasing the simulation length T 
is difficult: It becomes harder to preserve phase coherency, the 
outer boundary of a simulation is in causal contact for a larger 
fraction of the simulation, and existing codes would require 


* Besides its importance for GW astronomy, NR has also deepened the un¬ 
derstanding of general relativity in topics such as binary BH recoil [14, 15] 
and gravitational self force [16]. 


^ Several PN waveforms (or approximants) with different PN-truncation er¬ 
ror are available in the literature. These PN approximants can differ from 
each other during the last hundreds of cycles before merger. 
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TABLE I. Properties of the two NR simulations: The first block lists initial separation Dq, orbital frequency Hq, radial velocity do^ ADM 
energy Eadm angular momentum Tadm in units of total mass M = mj + m 2 . The middle block lists mass ratio mi Im 2 , eccentricity e, time 
duration T and number of orbits N until merger. The final block lists remnant mass Mf and spin 


many months or even years of wall-clock time. Therefore 
progress toward longer simulations has been sluggish, with T 
increasing by only about a factor of 2 to 3 during the last five 
years [25, 26, 28-30]. The duration T needed to close the gap 
depends on the binary parameters and the detector bandwidth. 
Here we start addressing the issue of the gap by focusing on 
the nonspinning case and high mass ratio, q = mi/m 2 = 7, for 
which the PN approximants can greatly differ [22, 31]. We 
present a new NR simulation that extends T by a factor of 20 
and reduces the initial frequency /ini by a factor of 3. With 
its comparatively high mass ratio, the new simulation probes 
an astrophysically relevant parameter regime for BH-BH and 
NS-BH binaries and for certain total masses covers the entire 
frequency band of advanced LIGO (aLIGO) and Virgo. We 
describe challenges involved in carrying out this new simula¬ 
tion, most notably an instability that causes the center of mass 
(CoM) of the binary to move, and we suggest improvements 
for future long simulations. We then compare the new simu¬ 
lation with existing analytical waveform models to assess the 
impact of waveform model errors on the detection rate of ad¬ 
vanced detectors. 

Numerical-relativity waveforms. We report on two new 
simulations of a nonspinning BH binary with mass ratio q — 
nil I m 2 = 7. The short simulation is of typical length: 20 
orbits, r = 4, lOOM. The long simulation, the main focus 
of this paper, is about 25 times longer. Both simulations 
are computed using the Spectral Einstein Code (SpEC) [32]. 
The short simulation uses established computational tech¬ 
niques [25]. The speed-up needed for the long simulation 
is the result of a series of code changes including task-based 
adaptive parallel load-balancing, live timing-based selection 
of the most efficient algorithm (when multiple implementa¬ 
tions of the same function are available), a modified mem¬ 
ory layout to allow more efficient calls to low-level numeri¬ 
cal packages and a more efficient implementation of the Gen¬ 
eralized Harmonic evolution equations. Figure 1 shows the 
new long waveform and Table I presents additional details 
about both simulations. Geometrized units G = c = I are 
used in Table I and throughout this paper. The top inset of 
Fig. 1 shows the spectra of the (2,2) spherical harmonic wave¬ 
form modes. The long simulation covers the entire design- 
aLIGO frequency range for nonspinning BH-BH binaries with 
M > 45 Mq, and covers the early-aLIGO frequency range for 
M > IIMq, including nonspinning NS-BH binaries^. In con- 
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FIG. 1. Overview of the new very long simulation. The main panel 
shows the (2,2) spherical-harmonic mode of the GW strain, with en¬ 
largements in the lower insets. The top inset shows the Fourier spec¬ 
tra of the new waveform in blue and the NR-NR hybrid waveform 
(used for comparisons with analytical models) in yellow, overlaid 
with noise power spectral densities of aLIGO at the early (dashed 
black) and design (solid black) sensitivity [4]. The waveforms in the 
inset are scaled to total mass M = 45.5Mq and luminosity distance 
Di ~ 1.06 Gpc. For comparison, an older q = 6 waveform [24] of 
representative length is shown in the main panel (offset vertically for 
clarity) and in the power-spectrum inset. 
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trust, the q = 6 simulation plotted in green, which is represen¬ 
tative of past simulations, starts at 3 times higher frequency, 
and covers a much smaller portion of the aLIGO bandwidth 
for a given M. Thus, we present here the first gravitational 
waveform covering the entire design-aLIGO frequency band 
for a nonspinning, compact-object binary at mass ratio q — 1 
with a total mass as low as M = 45.5M0. 

The short simulation is run at three different numerical 
resolutions, and the long one at four resolutions. The long 
simulation employs dynamical spectral adaptive mesh refine¬ 
ment [34], so measured quantities (like BH masses or wave¬ 
forms) do not always converge in a regular, predictable man¬ 
ner with increasing resolution, as is the case when each resolu¬ 
tion is defined by a fixed number of grid points. Furthermore, 


^ For mass ratio 7, in absence of spin, we expect no observable differences in the merger signal between a BH-BH and a NS-BH binary [33]. 
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FIG. 2. Top left: displacement of the CoM from the origin, for 
three long simulations with different outer boundary radii R. In each 
case, |ccomI increases exponentially. Top right: growth rate a of 
|ccomI as a function of R with a power-law fit. Lower panels: GW 
(2,2) mode of the short and long simulations. The long simulation 
still agrees very well with the short one at early times, but fails to 
produce an accurate merger waveform. 


failure to resolve initial transients caused by imperfect initial 
data also complicates convergence (see discussion in Sec. IIIB 
of Ref. [35]). Nevertheless, we find that differences in mea¬ 
sured quantities like waveforms, masses, and spins decrease 
rapidly with resolution, and Table I displays a conservative er¬ 
ror estimate obtained by taking the difference between the two 
highest resolutions. After the initial transients have decayed, 
we measure the mass ratio to be equal to 7 to within five sig¬ 
nificant digits, and the dimensionless spins to be < 10^®. The 
remnant mass and spin agree to within four significant digits 
between the short and long simulations. 

However, the long simulation encounters an unexpected 
problem. After a few 10,000M of evolution, the coordinate 
CoM begins to drift away from the origin. Define the CoM 
as (?CoM = c\m\/M + C 2 m 2 /M, where ci 2 are the coordinate 
centers of the apparent horizons of the two BHs. The top left 
panel of Fig. 2 shows that |ccom| increases exponentially with 
time, contradicting all known physical expectations. This drift 
is primarily a coordinate effect that only marginally affects 
most measurable quantities. For example, the linear momen¬ 
tum radiated to infinity, as computed from the waveform ob¬ 
tained by Cauchy-Characteristic extraction [36], is consistent 
with PN theory, and is too small to account for the motion 
of the CoM. To explore the drift in more detail, we repeat 
the long simulation with different values of R, the coordi¬ 
nate radius of the artificial outer boundary where we impose 
an outgoing-wave boundary condition. The top left panel of 
Fig. 2 shows the CoM drift for several values of R. We find 
that the exponential growth rate a behaves approximately like 
a as shown in the upper right panel of Fig. 2. For 

our standard choice of R = 864M, 1/a = 26,000M; this large 


timescale explains why the drift was not noticed in earlier, 
shorter simulations. 

We conjecture that the drift is caused by a coupling with 
the outer boundary. Such a coupling might arise through en¬ 
hanced reflections of the outgoing GW at the outer bound¬ 
ary. Our outgoing-wave boundary conditions [37] have the 
smallest reflection coefficient for spherical harmonic modes 
with small and a reflection coefficient of order unity when 
kR/i^ 1 [37], where k is the radial wavevector. With increas¬ 
ing |ccom| the emitted GW will have increasing high-^ content 
when decomposed on the outer boundary. 

From the bottom panel in Fig. 2, we see that the effects of 
the drift on the long waveform is confined to the last ^ 10 or¬ 
bits before merger. In the Fourier domain (see Fig. 1) one can 
clearly see the unusual behavior of the long waveform at high 
frequencies. Since these orbits are covered by the short wave¬ 
form, we hybridize the short waveform with the long one, thus 
replacing the final portion of the former. We adopt the hy¬ 
bridization method of Ref. [38]. We construct 9 NR-NR hy¬ 
brids by combining 3 versions of the long simulation and three 
versions of the short simulation. Each version may differ by 
the numerical resolution (which we label by “LevM^”, with 
larger integers indicating finer resolution), or by the degree 
of the polynomial used to extrapolate the waveform to infin¬ 
ity [36, 39] (which we label by “N,^”, where ^ is the poly¬ 
nomial degree.) In particular, we use (Lev3, N3), (Lev3, N2), 
and (Lev2, N3) for the long simulation, while we use (Lev5, 
N3), (Lev5, N2), and (Lev4, N3) for the short simulation. The 
fiducial NR-NR hybrid is built from the long (Lev3, N3) and 
the short (Lev5, N3) simulations; this pair of waveforms is 
blended over the interval t — fpeak G [—3252, —2252]M. In the 
top panel of Fig. I, we show in yellow the spectrum of the 
fiducial NR-NR hybrid. This spectrum behaves as expected 
close to merger, and is devoid of oscillations, just like the 
spectrum of the q = 6 simulation. Since we cannot estimate 
the impact of the coordinate drift on the phase error of the long 
waveform, we cannot make statements about the phase dis¬ 
agreement between the long waveform and analytical wave¬ 
form models. Nevertheless, we can compare the analytical 
models to the hybrid NR-NR waveform, and investigate how 
the results change when we vary the blending window where 
the hybridization is done. 

Comparison to analytical-relativity waveforms. Fig¬ 
ure 3 summarizes our comparisons between various an¬ 
alytical waveforms with the hybrid NR-NR wave¬ 

form h 22 - Shown is the unfaithfulness defined as 
# = 1 - max, 0,00 {h^f, h^f) /1 11 /1 1 ^ 22 ^ 11. Here t p and 

00 are the initial time and phase, ||/!|| = (h,h), and 

(huln) = 4Re/^“fii(/)/;5(/)/5„(/)d/, where 5„(/) is the 
zero-detuned, high-power noise power spectral density of 
aLIGO [40]. We consider the following analytical wave¬ 
form models from the LIGO Algorithm Library (LAL); the 
inspiral-only PN Taylor approximants [4 1 ] in the time domain 
(Taylor-Tl, T2, T4) and in the frequency domain (Taylor- 
F2), an inspiral EOB model (obtained from Ref. [20] by 
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FIG. 3. Unfaithfulness of the hybrid NR-NR waveform against several analytical waveform models. Left panel: inspiral-only comparisons. 
Right panel: IMR comparisons. Also shown in the left panel are SNR?jj,p/SNR|jjj (solid black line) and SNR^^^p/SNRj^jj (dot-dashed 

black line), and in the right panel SNR^jyjj^/SNR^jjj (solid black line). The blue area indicates the NR error. 


dropping any NR information, thus uncalibrated), the IMR 
EOBNR models that were obtained by calibrating the EOB 
model to NR simulations [20, 42, 43] (denoted in LAL as 
EOBNRv2, SEOBNRvl and SEOBNRv2), and the IMR phe¬ 
nomenological models that were built combining PN and 
NR results [17, 18] (denoted in LAL as PhenomB and Phe- 
nomC). All the time-domain IMR waveforms are tapered us¬ 
ing a Planck windowing function [44], both at the beginning 
and at the end. We generate the model waveforms start¬ 
ing from an initial GW frequency of = 0.01317. Eor 
inspiral-only models, we set /max = 0.01176/M, the mini¬ 
mum available final GW frequency among the time-domain 
Taylor models, a value close to the innermost-circular-orbit 
value in Schwarzschild spacetime (« 0.01083/M), whereas 
for the IMR comparisons /max = Quite interestingly, the 
inspiral-only comparisons give similar results when employ¬ 
ing directly the long simulation instead of the NR-NR hybrid. 

The blue shaded area in Pig. 3 represents the uncertainty 
in the NR waveforms, estimated by computing ,0' between 
the fiducial hybrid NR-NR waveform and the other 8 NR-NR 
hybrids. Because the inspiral-only and IMR curves are calcu¬ 
lated using different portions of the hybrid NR-NR waveform, 
the same model may have different values in the two panels 
for the same total mass. We vary the prescriptions used for 
the hybridization (namely, position and width of the blending 
window), and we hnd changes ^(10^^) in the unfaithfulness 
curves for low total masses. Thus, we consider our results ro¬ 
bust. If general relativity correctly describes the GW signals 
found in nature, then the unfaithfulness ^ plotted in Pig. 3 
yields a bound on the loss in detection rate due to modeling 
error. Eor sources uniformly distributed in space, the rela¬ 
tive loss in detection rate is ~ 3 (£/mm + <7e) (see Sec. VB in 
Ref [41]) where liMM is the minimal match of the template 
bank and ^7E=1 ~tnax^ (/j^^,/!^ 2 ^)/||/!^^||/||/! 2 |^|| is the inef¬ 
fectualness. Here A describes all the binary parameters, not 


just 00 and fo, and therefore c/e < Typically, c/mm = 3% 
in LIGO searches. Thus, to achieve < 10% loss in detection 
rate, it suffices that < 1% [41]. 

Quite remarkably, we hnd that the unfaithfulness of the un¬ 
calibrated inspiral EOB waveform is <0.1%, with a negligi¬ 
ble loss in detection rate due to modeling error. The agree¬ 
ment is of course better for the inspiral EOBNR waveforms 
(i.e., EOBNRv2, SEOBNRvl, SEOBNRv2) (# < 0.02%, 
left panel), and ^ < 0.2% for the IMR EOBNR waveforms 
(right panel). The closeness of all inspiral EOBNR waveforms 
strongly suggests that the different calibrations and variations 
in the dynamics and energy huxes of those EOBNR mod¬ 
els [20, 42, 43] do not impact the low-frequency part of the 
waveforms, but affect (in a minor way) only the last stages 
of the inspiral and the merger. The unfaithfulness of the 
time- and frequency-domain inspiral-only PN Taylor approx- 
imants varies between 0.1% and 10% depending on the bi¬ 
nary’s total mass and the PN approximant used In particu¬ 
lar, Taylor-T4, which has the best agreement with NR in the 
equal-mass case [30], has the largest disagreement with the 
new long q = 7 NR waveform. (The approximants Taylor-T2 
and Taylor-F2 are not displayed, but lie between Taylor-Tl 
and Taylor-T4). The PhenomB and C models were fitted to 
hybrids built with Taylor PN approximants and less accurate, 
short NR waveforms, which may in part explain the large dis¬ 
agreement we hnd. 

The new long NR waveform covers the entire design- 
aLIGO frequency band only for total mass M > 45.5Mq; 
for smaller M, the unfaithfulness calculations in Fig. 3 ne¬ 
glect the lowest frequency portion of the waveform visible 


The large unfaithfulness of some of the PN Taylor approximants is due to 
differences in the evolution of the frequency and its first time derivative 
during the late inspiral phase. 
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to aLIGO, down to ^ 10 Hz. To understand the signifi¬ 
cance of the missing GW cycles in the low-frequency por¬ 
tion of the bandwidth, we compute the signal-to-noise ratio 
(SNR) accumulated within the frequency range of our com¬ 
parisons (SNRinsp and SNRjmr for the left and right panel 
of Fig. 3, respectively), and compare it with the SNR accu¬ 
mulated over the entire inspiral (SNRfuii_insp) and the entire 
IMR (SNRfuii-iMR)- To cover the entire design-aLIGO band¬ 
width we use the calibrated FOB model of Ref. [20]. Suitable 
(squared) ratios of these quantities, which represent the frac¬ 
tion of total SNR that is accessible to our comparisons, are 
plotted in Fig. 3. These ratios are < 1 whenever GW cycles 
are missing in the range lOHz < / < or, in the case of 
SNR?,j,p/SNR|jj[_jjy[j^, also when the merger-ringdown signal 
is in band. We find that, even for total masses < 45.5M©, 
the unfaithfulness can still be a meaningful assessment of the 
quality of the analytical models, since a large fraction of SNR 
is accumulated. Because the merger-ringdown portion of the 
waveforms becomes increasingly important at higher masses, 
the inspiral-only comparisons at high mass cover only a small 
fraction of the entire SNR, as illustrated by the steep decline 
of SNR?„ 5 p/SNR|j[|_j|yjj^ in the left panel of Fig. 3. 

Conclusions. To detect and extract unique, astrophysical 
information from coalescing compact-object binaries, GW in¬ 
struments employ model waveforms built by combining an¬ 
alytical and numerical-relativity predictions [17-20]. Cur¬ 
rently, the main uncertainty of those waveform models is 
caused by the gap between the regimes of applicability of 
those methods. This uncertainty can be addressed and, even¬ 
tually, solved by running much longer NR simulations. In this 
work we start to tackle this issue by producing a BH-BH sim¬ 
ulation 20 times longer than previous simulations. Because 
of an unexpected drift of the CoM during the last 40 GW 
cycles, we construct the full NR inspiral, merger and ring- 
down waveform by hybridizing the long NR waveform with a 
new, short NR simulation. The hybrid NR-NR waveform cov¬ 
ers the entire band of advanced GW detectors for total mass 
> 45.5M0. Comparing to analytical waveform models, we 
find strong evidence that — at least for nonspinning binaries 
— the FOB formalism accurately describes the inspiral dy¬ 
namics in the so-far unexplored regime of 20 to 176 orbits be¬ 
fore merger, and combined with previous work [27] provides 
accurate waveforms beyond the limited range of calibration. 
Quite remarkably, the excellent agreement of FOBNR wave¬ 
forms holds also for uncalibrated inspiral FOB waveforms. 
PN-approximants have larger errors and more importantly the 
errors vary substantially depending on the specific PN approx- 
imant used. 
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